'''
Created on 2011-3-7

@author: summit
'''
import Image
import scipy
import scipy.ndimage
import scipy.io
import math
import matplotlib.pyplot as plt
import os

eps =  2.2204e-016
def bwdist(bw):
    '''
    bwdist is function of matlab image toolbox, here we implement it in python
    '''
    return scipy.ndimage.morphology.distance_transform_bf(1-bw)

def imgae2graydouble(img):
    '''
    Converts image to one channel (grayscale) double
    '''
    
    if img.dtype in [scipy.float32]: # image is a double
        if img.ndim == 3:
            img = scipy.array(Image.fromarray(img).convert("L"))
    else: # image is a int 
        if img.ndim == 3:
            img = scipy.array(Image.fromarray(img).convert("L"))
        img = scipy.cast[scipy.float32](img)
    return img

# whole matrix derivatives
def shiftD(M):
    return shiftR(M.transpose()).transpose()

def shiftL(M):
    # here , we should use reshape to get a column vector
    return scipy.concatenate((M[:,1:M.shape[1]], M[:, M.shape[1]-1].reshape(-1,1)), axis=1)

def shiftR(M):
    # here , we should use reshape to get a column vector
    return scipy.concatenate((M[:,0].reshape(-1,1), M[:,0:M.shape[1]-1]),axis=1)

def shiftU(M):
    return shiftL(M.transpose()).transpose()

def get_curvature(phi):  
    # notice we use float compute,so here you see minus 0.0
    dx=(shiftR(phi)-shiftL(phi))/2.0
    dy=(shiftU(phi)-shiftD(phi))/2.0
    dxplus = shiftR(phi)-phi-0.0
    dyplus = shiftU(phi)-phi-0.0
    dxminus = phi - shiftL(phi)-0.0
    dyminus = phi - shiftD(phi)-0.0
    
    dxplusy =(shiftU(shiftR(phi))-shiftU(shiftL(phi)))/2.0
    dyplusx =(shiftR(shiftU(phi))-shiftR(shiftD(phi)))/2.0
    dxminusy=(shiftD(shiftR(phi))-shiftD(shiftL(phi)))/2.0
    dyminusx=(shiftL(shiftU(phi))-shiftL(shiftD(phi)))/2.0
    
    nplusx = dxplus/scipy.sqrt(eps+(dxplus*dxplus )+(scipy.power((dyplusx+dy )/2.0,2)))
    nplusy = dyplus/scipy.sqrt(eps+(dyplus*dyplus )+(scipy.power((dxplusy+dx )/2.0,2)))
    nminusx= dxminus/scipy.sqrt(eps+(dxminus*dxminus)+(scipy.power((dyminusx+dy)/2.0,2)))
    nminusy= dyminus/scipy.sqrt(eps+(dyminus*dyminus)+(scipy.power((dxminusy+dx)/2.0,2)))
    
    
    return ((nplusx-nminusx)+(nplusy-nminusy)/2)

def showcontour(I, phi, i):
    
    plt.title('Evolution')
    #img = scipy.zeros((I.shape[0], I.shape[1], 3))
    #img[:,:,0] = img[:,:,1] = img[:,:,2] = I
    plt.imshow(I)
    plt.gray()
    plt.contour(phi, [0,0], colors='g')
    plt.title('%d Iterations'% i)
    

def simpleseg(I, init_mask, max_its, E, T, alpha):
    '''
    simple level set
    '''
    # ensures image is 2D double natrix
    I = imgae2graydouble(I)
    # create a signed distance map (SDF) from mask
    phi=bwdist(m)-bwdist(1-m)-0.5
    
    # main loop
    for its in range(1,max_its+1):
        D = E - abs(I-T)
        K = get_curvature(phi)
        F = -alpha*D + (1-alpha)*K
        
        dxplus = shiftR(phi)-phi-0.0
        dyplus = shiftU(phi)-phi-0.0
        dxminus = phi - shiftL(phi)-0.0
        dyminus = phi - shiftD(phi)-0.0
        
        gradphimax_x = scipy.sqrt(scipy.power(scipy.maximum(dxplus,0.0),2)-scipy.power(scipy.maximum(-dxminus,0.0),2))
        gradphimin_x = scipy.sqrt(scipy.power(scipy.minimum(dxplus,0.0),2)-scipy.power(scipy.minimum(-dxminus,0.0),2))
        gradphimax_y = scipy.sqrt(scipy.power(scipy.maximum(dyplus,0.0),2)-scipy.power(scipy.maximum(-dyminus,0.0),2))
        gradphimin_y = scipy.sqrt(scipy.power(scipy.minimum(dyplus,0.0),2)-scipy.power(scipy.minimum(-dyminus,0.0),2))
        
        gradphimax = scipy.sqrt(scipy.power(gradphimax_x,2)+scipy.power(gradphimax_y,2))
        gradphimin = scipy.sqrt(scipy.power(gradphimin_x,2)+scipy.power(gradphimin_y,2))
        
        gradphi = (F>0)*(gradphimax) + (F<0)*(gradphimin)
        
        # stability CFL
        dt = 0.5/abs(F*gradphi).max()
        
        # evolve the curve
        phi = phi + dt*F*gradphi
        
        # reinitialise distance function every 50 iterations
        if its % 50 == 0:
            phi = bwdist(phi<0)-bwdist(phi>0)
            
        # intermediate output
        if its % 20 == 0:
            showcontour(I, phi, its)
            
            filename = os.path.join("./r" ,str('%03d' % (its)) + '.bmp')
            plt.savefig(filename, dpi=100)
            plt.clf()
            
        
if __name__ == '__main__':
#    img = Image.open("liver.bmp")   # load the image
#    I = scipy.array(img)                # to the matrix
#    m = scipy.io.loadmat("mliver.mat")["m"]
    
#    plt.figure()
#    plt.subplot(2,2,1)
#    plt.imshow(I)
#    plt.title('Input Image')
#    plt.subplot(2,2,2)
#    plt.imshow(m)
#    plt.title('Initial Mask')
#    plt.show()
#    seg = simpleseg(I, m, 1000, 35, 170, 0.02)
#    plt.subplot(2,2,1)
#    plt.imshow(seg)
#    plt.title('Final Mask of phi<=0')
#    plt.show()

    import subprocess
    import glob
    filenames = glob.glob('./rv/*.png')
    filenames.sort(key=lambda filenames:int(filenames[5:-4]))
    filenames = [line+'\n' for line in filenames]
    print len(filenames)/10
    f = open('./list.txt', 'w')
    f.writelines(filenames)
    f.close()
    outname = "v"
    command = ('./mencoder',
           'mf://@list.txt',
           '-mf',
           'type=png:w=800:h=600:fps=%s'% str(len(filenames)/10),
           '-ovc',
           'lavc',
           '-lavcopts',
           'vcodec=mpeg4:mbd=2:trell:autoaspect',
           '-oac',
           'copy',
           '-o',
           './%s.avi' % outname)
    print "\n\nabout to execute:\n%s\n\n" % ' '.join(command)
    subprocess.check_call(command)

    
        
